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O ABSTRACT 

o ' 

Context. Understanding grain-surface processes is crucial to interpreting the chemistry in many regions of the interstellar medium. 
' However, accurate surface chemistry models are computationally expensive and are difficult to integrate with gas-phase simulations. 

Aims. A new modified-rate method for solving grain-surface chemical systems is presented. The purpose of the method is to trade a 
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small amount of accuracy, and certain excessive detail, for the ability to accurately model highly complex systems that can otherwise 
only be treated using the sometimes inadequate rate-equation approach. 

Methods. In contrast to previous rate-modification techniques, the functional form of the surface production rates was modified, and 
not simply the rate coefficient. This form is appropriate to the extreme "small-grain" limit, and can be verified using an analytical 
master-equation approach. Various further modifications were made to this basic form, to account for competition between processes, 
to improve estimates of surface occupation probabilities, and to allow a switch-over to the normal rate equations where these are 
applicable. 

Results. The new method was tested against a number of systems solved previously using master-equation and Monte Carlo tech- 
niques. It is found that even the simplest method is quite accurate, and a great improvement over rate equations. Further modifications 
allow the master-equation results to be reproduced exactly for the methanol-producing system, within computational accuracy. Small 
discrepancies arise when non-zero activation energies are assumed for the methanol system, which result from complex reaction- 
competition processes that cannot be resolved easily without using exact methods. Inaccuracies in computed abundances are never 
greater than a few tens of percent, and typically of the order of one percent, in the most complex systems tested. 
Conclusions. The new modified-rate approach presented here is robust to a range of grain-surface parameters, and accurately repro- 
duces the results of exact methods. Furthermore, it may be derived from basic approximations, making the behaviour of the system 
understandable in terms of physical processes rather than time-dependent probabilities or other more abstract quantities. The method 
is simple enough to be easily incorporated into a full gas-grain chemical code. Implementation of the method in simple networks, 
C*") . including hydrogen-only systems, is trivial, whilst the results are highly accurate. 
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ryQ 1. Introduction to be statistically large, making stochastic fluctuations unimpor- 

, tant. Such a "deterministic" treatment allows the chemistry to 

.. , Chemical models consisting of hundreds of species and thou- be desC ribed by a set of ordinary, first-order differential equa- 
t > • sands of reactions are frequently used to interpret the mor- tions (rate equations)i which ^ sorved numerically, allowing 
k> ; phologies and evolutionary histories of interstellar clouds and the time-dependent behaviour to be traced. 

, star-formation regions. Towards such ends, simulations deal- xhe simplicity of this rat e-equation approach naturally en- 
Ctf ■ ing solely with gas-phase chemistry have achieved much sue- courages its use in solving gra in-surface problems, especially 
' cess; however, the explicit consideration of grain-surface pro- - m coupled gas . p hase/grain-surface (i.e. gas-grain) models; how- 
cesses in chemical models is becoming increasingly important to ^ its accuracy in these applications is in some cases ques- 
our understanding of interstellar chemistry. From the formation t ionable. Grain-surface chemistry takes place on finite surfaces, 
of the most basic and abundant interstellar molecule, H 2 (e.g. where the populations of cert ain chemical species can become 
Gould & Salpeter 1963; Hollenbach & Salpeter 197 1 ; Duley & small of the order of j If surface reactions occur quickly 

Williams 1984), to some of the most chemically complex or- - m such a regim6i then stochastic effects may come int0 play; this 
game species observed in star-forming regions (e.g. Charnley et renders the rate . e q U ation treatment inaccurate, yielding produc- 
al. 1992, 1995; Cazaux et al. 2003; Horn et al. 2004; Garrod tion rates that are faster than physically possib le. 
& Herbst 2006; Garrod et al. 2008), the action of grain-surface These failures do not ^ in all regimes; | Ka tz et al. (1999)1 
chemical processes appears to be crucial. showed that hydrogen diffusion at low temperature takes place 

Unfortunately, the accurate coupling of gas-phase and grain- via thermal hopping, rather than fast quantum tunnelling. For 
surface processes in chemical models is difficult, due to the dif- grains of canon i ca i srze (O.ljum), at temperatures of -10 K, the 
ferent nature of the chemistry in each phase. Gas-phase chem- reac tion rates are typically no faster than the accretion or evap- 
istry may be accurately modelled by describing chemical abun- ora ti on rates of the reactants, putting the system into the deter- 
dances as averaged concentrations, or number densities. By as- m i n i st i c limit. It is in the consideration of grains that are smaller 
suming an arbitrarily large "cell" of gas, the absolute number of or hot ter than this that stochastic effects must be considered, 
particles of any species present in the system can be assumed Monte Carlo met hods have been developed (e.g. Tielens & 

Hagen 1982; Charnley et al. 1997; Charnley 1998, 2001) to 

Send offprint requests to: R. T. Garrod take account of stochastic behaviour. The most recent techniques 
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(Cuppen & Herbst 2005) trace the behaviour of individual atoms 
and molecules on a grain surface; however, even the simplest of 
such schemes is difficult to integrate with a rate equation-based 
gas-phase code (Chang et al. 2007). Alternatively, the master- 
equation method describes the system using the time derivatives 
of the probabilities of specific population states, which are eas- 
ily solved in tandem with gas-phase rate equations (Biham et 
al. 2001). Cut-offs are imposed, representing the maximum pop- 
ulations of surface species. Unfortunately, with large networks 
the number of possible population-state combinations becomes 
unmanageable, even with low cut-offs. 

Stantcheva et al. (2002) developed hybrid schemes that mix 
deterministic and stochastic methods, treating abundant, "deter- 
ministic" species such as CO using rate-equations, but treating 
"stochastic" species using a master-equation method that em- 
ploys low cut-offs. These schemes are successful for the small 
networks upon which they have been tested, but require prior 
knowledge from the exact methods in order to distinguish "de- 
terministic" from "stochastic" species. 

Green et al. (2001) made approximations within the master- 
equation method to give analytic expressions for production 
rates in some simple chemical systems; but the approximations 
are valid only in the stochastic regime, ruling out an extension to 
complex networks involving deterministic species such as CO. 

The "method of moments" recently developed by Barzel & 
Biham (2007a,b) involves the calculation of second moments of 
population, beginning from a master-equation standpoint. This 
allows the calculation of accurate production rates. However, in 
the deterministic regime, the method gives accurate production 
rates, but can produce inaccurate populations (Barzel & Biham, 
2007b; B. Barzel, private comm.). It is currently unclear how this 
would affect the behaviour of large chemical networks. Also, the 
scheme is untested against networks that include activation en- 
ergies, such as are necessary in the methanol-producing system. 

Attempts have been made previously (Caselli et al. 1998; 
Shalabiea et al. 1998; Stantcheva et al. 2001) to rectify the inac- 
curacies of the rate-equation method by modifying the rate co- 
efficients of grain-surface production rates according to a set of 
simple rules; Caselli et al. (2002) also devised a more complex 
modification scheme for reactions involving activation energy 
barriers. These modifications were largely empirical, and were 
not universally successful; convergence with exact techniques 
was achieved only for low temperatures, where atomic hydrogen 
is the only mobile reactant. The comparison of techniques con- 
ducted by Rae et al. (2003), using a standardised two-reactant 
system, found that the method was in general no more accurate 
than the standard rate-equation approach. 

Here is presented a new method of modification of the pro- 
duction rates, different from previous techniques. This method 
modifies not the rate coefficient, but the functional form of the 
production rate itself, adopting production rates derived for the 
extreme "small-grain" regime. The method switches smoothly to 
a rate-equation treatment where appropriate. The scheme retains 
its functional dependence on averaged population values, (AO, 
rather than analysing individual states, N, allowing it to be eas- 
ily implemented in standard gas-grain codes. The purpose of the 
new method is not primarily to produce a perfect match to the 
results of exact methods, but to improve accuracy in large chem- 
ical networks that cannot be modelled with exact techniques. 

Section 2 outlines the rates for surface processes. Basic mod- 
ified rates are formulated in Section 3. Sections 4 and 5 detail 
further modifications to the basic rates. Section 6 applies the 
method to reaction systems with activation energies. A discus- 
sion follows in Section 7, and conclusions in Section 8. 



2. Accretion, Evaporation and Reaction Rates 

In the systems to be analysed in this paper, three basic grain- 
surface processes are considered: accretion of gas-phase atoms 
onto the grains; evaporation of grain-surface species into the gas 
phase; and reaction between species on the grain surfaces. 

Here, the accretion rate is defined simply as an average flux 
of atoms/molecules per second. This obviates the need to trace 
any kind of gas-phase chemistry in the following tests. 

The evaporation rate of species A is governed by the lifetime 
against evaporation of one particle, t evap {A). Assuming a binding 
energy Ed(A) (in Kelvin), and a grain temperature Tj, the rate 
coefficient for evaporation is: 



k eV ap(A) = v • exp [-E D (A)/T d ] 



(1) 



where v is the vibrational frequency of the atom/molecule in its 



binding site, typically of the order of 10 s (e.g., Hasegawa et 
al. 1992). Hence the total evaporation rate (in s~') is: 



RevapfA) = k emp (A) ■ (N(A)) 



(2) 



where (N(A)) is the average surface population of species A. 

It is assumed that surface reactions take place via the 
Langmuir-Hinshelwood mechanism: reaction occurs when two 
reactants diffusing over the grain surface meet in the same bind- 
ing site. The rate at which reaction occurs is: 



_ Rhop(A) + Rhop(B) 

K-AB - ^ Kab 



(3) 



where Ri lop (A) and Ri, op (B) are the average rates of thermal hop- 
ping between binding sites for each reactant, S is the number of 
binding sites on the grain, and kab is an efficiency factor that ac- 
counts for any activation energy required to react. If diffusion is 
assumed to occur via thermal hopping, the rate is simply 



R ho p(A) = vexp[-E b (A)/T d ] 



(4) 



where E/, is the strength of the barrier (in Kelvin) between bind- 
ing sites. If, in the case of atomic hydrogen, movement between 
sites is assumed to occur through quantum tunnelling, different 
expressions pertain; see Hasegawa et al. (1992)[ 



Where an activation energy is required for reaction to occur, 
the efficiency, kab> is typically described by a Boltzmann factor, 
or an expression for quantum-tunnelling efficiency, whichever 
is the more efficient. Alternatives to these simple expressions 
have been put forward that take account of competition between 
reaction and diffusion out of the binding site (Awad et al. 2005, 
Chang et al. 2007); however, for the purposes of comparison to 
older results, these more recent developments are ignored. 

Using rate equations to solve a grain-surface chemical net- 
work, the production rate (in s~') of the products of the reaction 
between species A and B is defined as: 



RprodiAB) = k AB ■ (N(A)) ■ (N(B)) 

or, in the case of homogeneous reactants, 



RproAAA) = J± ■ (N(A)f. 



(5) 



(6) 



As stated above, the accuracy of this definition breaks down 
when stochastic effects become important. 

The difference between the production rate, R pro d(AB), and 
the reaction rate, k AB , should be noted. The former is the subject 
of the modifications proposed in this paper; the latter is strictly 
defined above in equation (3), and is not adjusted in any way. 
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Fig. 1. Population of H atoms and production rate of H2 on 
grain surface as a function of S , the number of sites per grain, 
adopting low R acc (H). Red lines are the master-equation results 
of |Barzel & Biham (2007b) with crosses indicating the data 
points. Solid black lines indicate results of the simple modifi- 
cation scheme of Section 3.1. Dashed lines are the standard rate- 
equation solutions; dotted lines indicate the solutions using only 
equations (11) and (12). 




10 3 10 4 10 5 10 6 



Fig. 2. Population of H atoms and production rate of H2 on 
grain surface as a function of S , the number of sites per grain, 
adopting high R acc {H). Red lines are the master-equation re- 
sults of Barzel & Biham (2007b) , with crosses indicating the 
data points. Solid black lines indicate results of the simple mod- 
ification scheme of Section 3.1. Dashed lines are the standard 
rate-equation solutions; dotted lines indicate the solutions using 
only equations (11) and (12). 



2.1. Rate-equation inaccuracies 

The expectation value, (N(i)), of the grain-surface population of 
species i may be defined simply as: 

(N(i)) = J^N-P N (i) (7) 

N=0 

where fV(0 is the probability of finding N atoms/molecules of 
species i on the grain, at any arbitrary moment. The quantity 
(iV(0) may be considered to be the average population of species 
i over a large ensemble of interstellar dust grains, as suggested 
by |Lipshtat et al. (2"0 04); but it may also be viewed as the av- 
erage population on an individual grain over a statistically long 
time period. Such an understanding is valid so long as the macro- 
scopic behaviour of the system is in a quasi-steady state in rela- 
tion to the microscopic processes that occur on grains. 

As explained by Biham et al. (2001) and Lipshtat et al. 
(2004), the inaccuracies of rate equations result from the use of 
the multiplication products (N(A)) ■ (N(B)) and ±(JV(A)> 2 in the 
production rates of equations (5) and (6). These quantities are 
used to represent the average number of unique pairs of reac- 
tants present on the grain. They are thus only approximations 
to the true average number of unique pairs, (N(A) ■ N(B)) and 
l(N(A) ■ [N(A) - 1]). The approximations are valid when the 
surface populations of reactants are large, but may break down 
when those values become small. Hence, the equalities, 

(N(A)) ■ (N(B)) = (N(A) ■ N(B)) (8) 



(N(A)f = (N(A) ■ [N(A) — 1]) = (N 2 (A)) - (N(A)) (9) 

are expressions of the validity of the "deterministic" rates shown 
in equations (5) and (6). However, the condition that (N(A)) or 
(N(B)) be small (i.e., of the order of 1 or less) is necessary, 
but not sufficient, to invalidate equations (8) and (9). This may 
clearly be seen from the results of Barzel & Biham (2007b), in 
which the deterministic rates are sometimes seen to be accurate 
even when the populations of all reactants are comfortably less 
than 1 . In the case of equation (8), the equality is valid so long 
as the probabilities, Pn(A) and Pn(B), of individual population 
states, N(A) and N{B), remain uncorrected. If average popula- 
tions are small and reaction rates, given by equation (3), are fast, 
then these probabilities may become anf/-correlated, invalidat- 
ing equation (8). In this case, the probability of finding one par- 
ticle of A and one particle of B on the grain surface at the same 
time is low, because they quickly react to produce AB. Following 
this argument, the deterministic rates, equations (5) and (6), must 
represent an absolute upper limit to the production rates, in all 
regimes. 

3. Basic modifications 

Here, rather than calculate (N(A) ■ N(B)) or (N 2 (A)), the grain- 
surface production rates are evaluated by an altogether different 
approach. Firstly, a basic form is constructed for the production 
rates, with which the standard reaction rates of equations (5) and 
(6) may be replaced. In order to retain the same matrix-inversion 
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Table 1. The surface reactions used in each system. Activation energies employed by Stantcheva et al. 2002 are also indicated. 



Reaction 


Barzel & Biham (2007b) 


Stantcheva et al. 


(2002) 


T T j_ T T A . /"IT T /M T j_ 

H 2 system H 2 system CH 3 OH system 


CH3OH system 


E A (K) 


H + H -> H 2 


• • • 


• 




H + O -> OH 


• • 


• 




H + OH -> H 2 


• • 


• 




O + o -» o 2 


• • 


• 




H + CO -> HCO 




• 


2000 


H + HCO -» H 2 CO 




• 




H + H 2 CO -> CH 3 




• 


2000 


H + CH3O -> CH3OH 




• 




O + CO -> co 2 




• 


1000 


O + HCO — > C0 2 + H 




* 





techniques with which the gas-phase equations are integrated, it 
is essential to express these rates in terms of the averaged popu- 
lations, (N(i)), of each species ;. 

Consider a system of two reactive species, A and B, accreting 
onto a dust grain at rates R aC c(A) and R acc (B), and with evapora- 
tion rates k evap (A) and k emp (B). A single reaction, A + B — » AB, 
is allowed to occur, at a rate k AB , as defined in equation (3). 

The simplest stochastic case is that in which each reacting 
species has an average abundance (N(i)) « 1. Here, the prob- 
ability of population states N(i) > 2 is very small, so the prob- 
ability of finding 1 (or more) of reactant i on the grain at any 
moment may be approximated as: 

P(0 - (N(i)). (10) 

It is assumed that reaction between A and B is fast; i.e., 
k A B » Racc(A),R acc (B),k emp (A),k evap (B), so A and B will re- 
act as soon as both are present. The reaction rate is therefore 
unimportant in the calculation of the overall production rate of 
species AB (thus the rate-equation method fails). It is the ac- 
cretion rates that determine the production rate; the system is 
in what is commonly referred to as the "accretion limit". The 
rate of production is simply equal to the rate of accretion of one 
species times the probability that the other is present, and vice 
versa. Using equation (10) this gives: 

R mod (AB) = R acc (B) ■ (N(A)) + R acc (A) ■ (N(B)). (11) 

If A and B are atoms/molecules of the same species (e.g. 
atomic H), which may react together, the final result is: 

R mod (AA)=R acc (A)-{N{A)). (12) 

As long as k AB » k evap (A), k evap (B), there will be no com- 
petition from other processes; hence, expressions (11) and (12) 
are accurate even if the reaction rate is comparable to either ac- 
cretion rate. Implicit in the formulation above is the assumption 
that only the two species crucial to the reaction are present on 
the grain; no competition with other reactions is considered. 

Equations (11) and (12) may also be easily arrived at by 
an analytical master-equation approach. Lipshtat et al. (2004) 
dub this result the "small grain" approximation, as low values of 
(N(A)) and (N(B)) and fast reaction rates, k AB , are achieved as 
the grain radius asymptotically approaches zero. 



Expressions (11) and (12) provide the basic form of the pro- 
duction rates with which it is proposed to replace the old rate- 
equation expressions, given the correct conditions. Equations (5) 
and (6) may be regarded as representative of the extreme situa- 
tion (N(A)), (N(B)) » 1, and equations (11) and (12) as repre- 
sentative of the other extreme, (N(A)), (N(B)} « 1. 

This simple formulation may be compared to that proposed 
by Caselli et al. (1998) and explored in a number of subse- 
quent papers. They proposed, in their "corrected" formulation 
(Stantcheva et al. 2001), to replace the standard hydrogenation 
reaction rates, k^x, i-e. the production rate coefficients, with the 
faster of either the accretion rate or evaporation rate of atomic 
hydrogen. In the case where k evap (H) > R acc (H) (the so-called 
"evaporation limit"), the steady-state abundance of atomic hy- 
drogen may be evaluated as (iVTH)) ^ R acc (H)/k emp (H), yielding 
a modified production rate of R acc (H) ■ (N(X)). For the reaction H 
+ H — > H2, this form is equal to equation (12). For other reactions 
(with X^H), (N(H)) is expected to be very small, so this same 
form is also approximately equal to equation (11). Hence, the 
modifications for hydrogenation reactions proposed by Caselli 
et al. (1998) and Stantcheva et al. (2001) for the "evaporation 
limit" case may be seen to be a special case of equations (11) and 
(12). Those modifications remain accurate so long as the "small 
grain" limit is maintained, and evaporation is still the dominant 
removal mechanism for atomic hydrogen. 

3.1. Rate replacement and restrictions 

To determine whether the production rates of equations (5) and 
(6) should be replaced with equations (11) and (12), the validity 
of equations (8) and (9) must be assessed. If (N(A)} » 1 and 
(N(B)) » 1 then the equalities, and the deterministic rates, are 
valid. In the case that, for example, (N(A)) » 1 but (N(B)) « 
1, equation (8) is still valid, and rate equations may be used. 
(An assumption to this effect was made by Stantcheva et al. 
2002, in treating species such as CO using deterministic rates). 
In such cases, reactions involving species B would have only a 
small effect on the population state of species A, making any 
(anti-)correlation very weak. In the case where both (N(A)) « 
1 and (N(B)) « 1 the rate substitutions should be made. 

To apply these conditions, species with (N(i)) > 1 are 
deemed to be in the (N(i)) » 1 regime, whilst species with 
(N(t)) < 1 are deemed to be in the {N(i)) « 1 regime. 
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Table 2. Binding energies, diffusion barriers, fluxes, and dust temperatures employed by Barzel & Biham (2007a,b) 
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Species 


H 2 system 




H 2 system 




CH 3 OH system 




(T d = 10 K) 




(T d = 15 K) 




(T d = 15 K) 




t D (K) t b (K) Racc/il (S ) 


fifl (K) 


E b (K) 


RaccIS (s- 1 ) 


Co (K) 




RaccIS (S" 1 ) 


H 


371 255 1.00 X 10 2.75 X 10 


603 


511 


5.0 x 10-'° 


603 


511 


5.0 x 10- 10 


O 




627 


545 


1.0 x 10-" 


627 


545 


1.0 x 10-" 


OH 




627 


545 




627 


545 




CO 










638 


580 


1.0 x 10- 10 


HCO 










673 


603 




H 2 CO 










685 


615 




CH 3 










719 


638 





S equals the number of binding sites on the grain 



Modifications are therefore made only when (N(A)), {N(B)) < 1 . 
This is, of course, a gross simplification; however, the tests to 
follow demonstrate that any resultant inaccuracies are small. 

As a further restriction on the new rates, it is asserted that 
under no circumstances may the modified production rate ex- 
ceed the standard rate-equation value, following the argument of 
Section 2. 1 . It is therefore required that: 



Rmod(AB) < k AB ■ (N(A)) ■ (MB)}. 



(13) 



Thus, if a modified production rate exceeds the deterministic 
rate, the deterministic rate is used. 



3.2. The hydrogen system 

This basic formulation is tested against the simplest of grain- 
surface systems, in which atomic hydrogen is the only reac- 
tive species. Barzel & Biham (2007b) investigated this sys- 
tem for various grain sizes, at a temperature of 10 K. They 
conducted rate-equation, master-equation, and moment-equation 
simulations, run to steady state in the population of atomic hy- 
drogen. Calculations were made using either a high or a low 
flux of accreting H-atoms; R aC c(H) = 2.75 x 10~ 8 S s -1 or 
Racc(H) = 1 x 10 _11 5 s , where S is the number of surface 
binding sites. Details of this system are given in Tables 1 & 2. 
|Barzel & Biham (200 7b) detail all other relevant values. 

Solving for H populations and H2 production rates requires 
the solution of the equation: 



d(N(H)) 
dt 



= R acc (U) - k evap (U) ■ <JV(H)> - R prod (H 2 ) = 



where R pro d(H2) is defined according to the stipulations of 
Section 3.1. Obtaining an analytical solution is trivial. Figures 
1 and 2 show grain-surface atomic hydrogen populations and 
H2 production rates calculated in this way for various values 
of S, using low and high H-fluxes, respectively. Solid black 
lines show the results when modified rates are employed accord- 
ing to the stipulations of the modification scheme of Section 
3.1. Dashed lines represent standard rate-equation results; dot- 
ted lines represent results obtained purely from equations (11) 
and (12). Indicated in red are the master-equation results of 
|Barzel & Biham (2007b)| with crosses marking the individual 
data points. The master-equation results may be regarded as a 
true representation of the system. 



Rates and populations calculated using standard rate equa- 
tions rise linearly with increasing S . At large grain sizes these 
are the exact solutions, but become inaccurate for small grains. 
Similarly, populations calculated purely with equations (11) and 
(12) may diverge from the rate equations for large grain sizes. 

The simple modification scheme demonstrates near-perfect 
agreement with the master-equation populations and produc- 
tion rates calculated by Barzel & Biham, at both the large- 
and small-grain extremes of the system, for each H-flux value. 
Approaching the stochastic-deterministic threshold, where de- 
terministic rates become accurate, results vary marginally from 
the exact master-equation values; the match is otherwise perfect. 

4. Continuous modification schemes 

Whilst modification using equations (10) - (13), along with the 
stipulations of Section 3.1, is accurate in the case of simple, 
highly prescribed systems, the purpose of this work is to fash- 
ion a modification scheme that is more universally applicable 
to large gas-grain chemical networks. Such a scheme must pro- 
vide rate continuity over the stochastic-deterministic threshold 
at (N(i)) = 1. Any functional form must also be integrable using 
standard differential equation-solving techniques (i.e. the Gear 
algorithm). Below, further modifications to, and restrictions on, 
the basic equations (11) and (12) are formulated. 

4.1. Threshold continuity 

The switch-over between the modified rate and the standard 
rate may produce a discontinuity in the rates at (N(i)) = 1. 
Dependent on the relative rates of all the processes involved, this 
may present an impediment to accurate calculations. Therefore, 
a simple empirical function, /, is introduced to make the tran- 
sition smoother whilst preserving a fast switch-over. Under this 
scheme, production rates are always modified according to: 

Rtot = Sab ■ RmodiAB) + (1 - f AB ) ■ k AB ■ (N{A)) ■ (N(B)) (14) 

where: 



(N(A))<l,(N(B))<l 
(N{A))>l,(N(B))<l 
{N(A))<1,(N(B))>1 
(N(A))>l,(N(B))>l 



/ab = 1 

f AB = 1/<JV(A)> 
fAB = V(N(B)) 
f AB = l/[(N(A))-(N(B))l 



(15) 
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Fig. 3. Populations and production rates for the H2O system, 
assuming no reaction competition. Red lines are the master- 
equation results of Barzel & Biham (2007b), with crosses indi- 
cating the data points. Solid black lines are the results of the 
modification scheme of Section 4. 1 . Dashed lines indicate the 
standard rate-equation results. 



Fig. 4. Populations and production rates for the CH3OH sys- 
tem, assuming no reaction competition. Red lines are the master- 
equation results of Barzel & Biham (2007b) with crosses indi- 
cating the data points. Solid black lines are the results of the 
modification scheme of Section 4.1. Dashed lines indicate the 
standard rate-equation results. 



This has the effect that when (N(A)),(N(B)) < 1, the rate 
is always "stochastic", whilst quickly tending towards the de- 
terministic rate as either (N(A)) or (N(B)) rises above unity. 
For example, if species A and B both attain abundances of 10 
atoms/molecules per grain, the deterministic contribution to the 
total production rate will be 99% of the normal deterministic 
rate. Expressions (14) & (15) therefore allow the modified rate 
contribution to replace a fraction of the total deterministic rate 
that corresponds to the reaction of 1 atom/molecule of species 
A with 1 atom/molecule of species B. Equation (13) is also ap- 
plied, meaning that if R mod (AB) > k AB ■ (N(A)) ■ (N(B)), the total 
rate is equal to the unmodified deterministic rate. 

4.2. The water system and the methanol system 

The continuous rate-modification scheme outlined above is used 
to examine the water and methanol systems, at 15 K, as investi- 
gated by Barzel & Biham (2007b) The former scheme includes 
reactions between surface species H, O, and OH, resulting in 
H2, O2 and H2O production. The latter scheme also includes re- 
actions with CO, HCO, H 2 CO, and CH3O, leading to the pro- 
duction of methanol, CH3OH, and carbon dioxide, CO2. The re- 
actions, fluxes, and binding energies are indicated in Tables 1 & 
2. It should be noted that the methanol system investigated by 
Barzel & Biham does not include any activation-energy barri- 



ers. In fact, activation energies substantially complicate the be- 
haviour of the methanol system; see Section 6. 

Figures 3 & 4 show population sizes and production rates of 
key species, for the water and methanol systems, respectively. 
Populations obtained from the new method (Section 4.1) are 
very well matched to the master-equation results of Barzel & 
Biham. Production rates are accurate at high and very low val- 
ues of S , but values near the stochastic-deterministic threshold 
are more obviously inaccurate, whilst obeying the correct trend. 
However, the new method is a very good first approximation. 

What are the underlying physical reasons for the disagree- 
ments? The production rates obtained from the new method rise 
to the rate equation values at lower values of S than the master- 
equation results. This occurs before any reactants approach a 
population of 1, as the rate limit of equation (13) is reached be- 
fore this. Hence, the empirical function that is used to switch 
over at the (N(i)) = 1 threshold is not the cause. 

In fact, the modified rates are too fast because competi- 
tion between surface processes has not been considered. At low 
values of S, reactions are extremely fast, due to the fast re- 
action rates of all reactants; see equation (3). As S increases, 
reaction becomes slower, and the possibility arises that one or 
other reactant may evaporate before the two meet in the same 
binding site and react. For hydrogen atoms in the water and 
methanol systems, using the binding energies shown in Table 



R. T. Garrod: A new modified-rate approach for gas-grain chemical simulations 



7 




10 2 10 3 10" 10 5 10 2 10 3 10" 10 5 

s s 



Fig. 5. Production rates for the H2O system (left panel) and CH3OH system (right panel), using the simple competition scheme of 
Section 5. Red lines are the master-equation results of Barzel & Biham (2007b), with crosses indicating the data points. Solid black 
lines show the results of the modification scheme of Section 5. Dashed lines indicate the standard rate-equation results. The apparent 
small deviations between master-equation and modified-rate methods derive largely from differences in S -sampling. 



2, knx - k evap (H) when S 460 sites per grain. Discrepancies 
in hydrogen-related production rates are greatest close to this 
point, as may be seen in figures 3 & 4. To accurately evaluate 
the production rates, competition processes must be considered. 

5. Evaporation-Reaction Competition 

The production rates expressed in equations (11) & (12) are sim- 
ply the rates at which one or other reacting species accretes onto 
the grain when the other is present on the surface. This assumes 
100% efficiency in the reaction that follows the accretion event. 

To take account of competition in the modified rates, an effi- 
ciency factor, Cab, is employed, such that: 

RmodiAB) = C AB ■ (Racc(B) ■ (N(A)) + R acc (A) ■ (N(B))) (16) 

Cab = k AB /[k A B + k evap (A) + k evap (B)]. (17) 

The efficiency, Cab, represents the probability that, of three pos- 
sible processes - reaction of A and B, evaporation of A, and 
evaporation of B - it is reaction that occurs. This assumes only 
two reacting atoms/molecules to be present at any one time; an 
assumption made by many other authors in considering such sys- 
tems (e.g. Allen & Robinson 1977, Green et al. 2001, Barzel & 
Biham 2007). Hence, reactions with other (stochastic) species 
may be ignored, for the purposes of competition. 

It should be noted that, when valid, the rate equations natu- 
rally take account of any form of competition that may arise. 

Applying the efficiency factor to R mo d(AB) for all reactions 
in the ITO and CH3OH systems gives the results shown in fig- 
ure 5. Production rates now show precisely the behaviour ex- 
pected, and provide an excellent match to the master-equation 
results, within computational accuracy. (Population sizes do not 
vary noticeably from the values shown in figures 3 & 4, and 
are hence omitted). Some small discrepancies seem to appear in 
the production rates, particularly for the CH3OH system; how- 
ever, these derive largely from greater S -sampling in the new 
modified-rate models than in the master-equation models. 

5.1. Anterior Competition and Formation Rates 

In the simple hydrogen system, the only source of hydrogen 
atoms on the grains is direct accretion from the gas phase. 



However, in the methanol system, hydrogen atoms may be pro- 
duced by chemical reaction. Furthermore, important radicals like 
HCO are formed solely on the grains, with no contribution from 
accretion. Thus, for the reaction H + HCO — » H2CO, equation 
(12) contains the single term fl flce (H) • (TV(HCO)). But, should 
HCO formation by reaction also be thought of as an "accre- 
tion" process, for the purposes of equation (12), leading to a term 
R f orm (HCO) • <tf(H)>7 

The inclusion of the term R / ofm (HCO) • {N(H)) would as- 
sume an H-atom to be present on the grain whilst HCO itself 
were formed. However, since HCO is formed by the reaction H 
+ CO — > HCO, this would mean the presence of two H-atoms on 
the grain at the same time, prior to HCO formation. Rather than 
forming HCO, the most likely outcome is that the two hydrogen 
atoms should react together instead. To treat this situation ac- 
curately, competition between reactions prior to the reaction of 
interest (i.e. H + HCO — » H2CO) should be taken into account. 
Such competition is labelled here as "anterior competition". 

In fact, if the two-particle approximation (e.g. Section 5) is 
valid, the problem of anterior competition may be almost en- 
tirely avoided. The two-particle assumption precludes the possi- 
bility that one reactant may be formed on the grains in the pres- 
ence of the other. This allows the formation processes resulting 
from stochastic reactions to be omitted from equations (11) and 
(12), or equation (16), and anterior competition ignored. 

However, in the case where, for example, CO is present in 
great abundance on the grains, the reactions of CO would be 
deterministic in nature, and CO should be considered present at 
all times. Thus, the deterministic formation of HCO by surface 
reactions should be included in the "accretion" rates utilised in 
equations (11) and (12), or equation (16). 

In the models that follow, R acc is substituted with Rf orm . 

R form (A)=R acc (A)+ £ (I - f u ) ■ kir (N(i)) ■ (N(j)} (18) 

all ij pairs 

for all reaction pairs ij that form species A, where kij is the re- 
action rate, and fij is defined in Section 4. 1 . As the deterministic 
part of the production rate in equation (14) is unaffected by the 
stochastic part, there is no need to iterate the calculations. 

For consistency, equation (17) should include, in the denom- 
inator, terms for the reaction of particles A and B with other 
species, such as CO, that can build up significant abundances 
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Table 3. Binding energies, diffusion barriers, and fluxes employed by Stantcheva et al. (2002) 



Species 


E D (K) 


E b (K) 




Rare (S _1 ) 




Low fin 


Intermediate n# 


High n H 


H 


350 


100 


1.67 x 10~ 5 


1.67 x 10~ 5 


1.60 x 10~ 5 





800 


240 


3.26 x 10~ 7 


2.72 x 10~ 6 


2.53 x 10~ 5 


OH 


1260 


378 


- 


- 


- 


H 2 


450 


135 


- 


- 


- 


o 2 


1210 


363 


- 


- 


- 


H 2 


1860 


558 








CO 


1210 


363 


2.05 x 10~ 7 


2.05 x 10~ 6 


2.05 x 10" 5 


HCO 


1510 


453 








H 2 CO 


1760 


528 








CH 3 


2170 


651 








CH3OH 


2060 


618 








co 2 


2500 


750 









on the grains. If CO, or any other species, has (N(i)} » 1 then 
there is a vanishing probability of not finding such a reactant on 
the grains when two other particles are also present. Reactions 
with that species must be allowed to compete with the reaction 
between A and B, if such reactions exist. 

In this model, (N(i)) = lis used as the stochastic- 
deterministic threshold. If (N(i)) > 1 then terms of the form 
(k-,j ■ [{N{i)) - 1]) should be inserted into equation (17) for the 
reactions of any species j with which species i may react. This 
means that only the deterministic part of the reaction rate is in- 
volved in the competition. This treatment is in keeping with the 
formalism of equations (14) and (15) and Section 4.1. 

In fact, this final precaution actually has no significant effect 
on the systems modelled in this paper; the activation energy bar- 
riers used later, in Section 6, are too high to make reaction with 
CO or H2CO sufficiently competitive to affect the reactions of 
atomic hydrogen or oxygen. However, for adoption in a gener- 
alised system, this eventuality is easily treated. 

6. Systems with activation energies 

Activation energy barriers are typically assumed to mediate key 
reactions in the methanol system. However, these barriers were 
not considered in the study of Barzel & Biham (2007b) In order 
to test the new modified-rate method against a system with acti- 
vation energies, the methanol system of Stantcheva et al. (2002) 
is employed; see Tables 1 and 3. 

|Stantcheva et al. (2002)| used a hybrid master-equation/rate- 
equation method to reproduce Monte Carlo results. The former 
approach assumed that the very reactive species H, O, OH, HCO 
and CH3O behave stochastically, but that species with reactions 
requiring activation energies (CO and H2CO), which typically 
achieve large populations, could be treated using the standard 
rate equations. The cut-offs for stochastic species were set to 
values 1, 2, or 3. Three basic scenarios were considered, charac- 
terised by the accretion fluxes of H, O, and CO atoms/molecules, 
mimicking conditions in interstellar clouds of "low", "interme- 
diate", and "high" density. The models were run to 1000 yr of 
evolution, assuming 10 6 binding sites per grain and a grain tem- 
perature of 10 K. Importantly, Stantcheva et al. allowed diffu- 



sion of hydrogen atoms around the grain surface to occur via 
quantum tunnelling between binding sites, rendering the rates of 
hydrogen reaction substantially faster than would be the case as- 
suming only thermal hopping. |Katz et al. (19 99) suggested that 
such quantum-tunnelling effects do not govern surface diffusion 
rates. However, the purpose of the comparison is, in this case, to 
test the method against a reliable standard, rather than to provide 
an accurate reproduction of grain-surface chemistry in interstel- 
lar clouds. More rigorous application of the new modification 
methods to interstellar cloud conditions will follow in future. 

Of the two methods utilised by Stantcheva et al., the Monte 
Carlo results should be considered the most reliable, although 
the differences are typically not large. However, the Monte Carlo 
simulations use only integer values, so it is not possible to obtain 
accurate estimates of average population sizes when those values 
are close to unity. Comparison with species of very low abun- 
dance relies on the results of the master-equation/rate-equation 
hybrid method alone. For convenience, this method will from 
now on be refered to simply as the "master equation" method, 
with the understanding that it is in fact a hybrid scheme. 

Table 4 shows two sets of results from Stantcheva et al.; the 
column headed "Master Equation 222 11" indicates results from 
their model using cut-offs for stochastic species H, O, OH, HCO, 
CH3O of 2, 2, 2, 1, 1, respectively. Columns in tables 5 and 6 are 
labelled similarly. Stantcheva et al. also ran models with lower 
cut-offs; the results of the models with the highest cut-offs are 
quoted, on the assumption that these are the most accurate. 

Tables 4-6 also show rate-equation solutions for each 
regime. Whilst a few rate-equation results are within an order 
of magnitude of the exact value, many populations vary wildly 
from the Monte Carlo and master-equation results, including 
important species such as formaldehyde, methanol and CO2. 
In fact, the rate-equation results shown in Tables 4-6 differ 
from the values quoted by Stantc heva et al. (2002)| which actu- 
ally correspond to activation energies of 2500 K, for the H + 
CO and H + H 2 CO reactions, rather than the stated 2000 K. 
Nevertheless, it is clear that the rate equations produce a very 
poor match to the exact solutions of each system. 

In Table 4, the column headed "Method A" shows low- 
density results using the simple competition scheme of Section 
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Table 4. Modified rate results for low density, using input values from Stantcheva et al. (2002). Values in boldface show agreement 
within 10% of the Monte Carlo or master-equation values. 



Species 


Stantcheva et al. 


Rate Eq. 




Modified Rates 




Monte Carlo 


Master Eq. 22211 




Method A 


Method B 


Method C 


Method F) 


"Simple" 

Cnmnpti ti on 


"Full" 
Competition 


Method B + 
Poisson Prob. 


Method C + 

1-T- accretion 


H 


1 


7.96(-3) 


1.22(-5) 


7.97(-3) 


7.97(-3) 


7.97(-3) 


7.96(-3) 


O 





1.90(-2) 


5.22(-7) 


1.85(-2) 


1.85(-2) 


1.86(-2) 


1.88(-2) 


OH 





1.93(-2) 


5.22(-7) 


1.86(-2) 


1.86(-2) 


1.88(-2) 


1.90(-2) 


H 2 


2 


1.94(+0) 


l.ll(+2) 


1.94(+0) 


1.94(+0) 


1.93(+0) 


1.93(+0) 




1.70(+2) 


1.62(+2) 




1 SQYj-? , | 

1 .07^TiJ 


1.89(+2) 


1.89(+2) 




H 2 


9.90(+3) 


9.86(+3) 


1.03(+4) 


9.76(+3) 


9.76(+3) 


9.76(+3) 


9.85(+3) 


CO 





2.81(-2) 


1.85(+1) 


3.76(-2) 


2.82(-2) 


2.82(-2) 


2.82(-2) 


HCO 





1.22(-2) 


3.29(-7) 


1.21(-2) 


1.21(-2) 


1.22(-2) 


1.22(-2) 


H 2 CO 





2.83(-2) 


1.89(+1) 


3.74(-2) 


2.82(-2) 


2.82(-2) 


2.84(-2) 


H 3 CO 





1.23(-2) 


3.29(-7) 


1.21(-2) 


1.21(-2) 


1.22(-2) 


1.22(-2) 


CH3OH 


6.40(+3) 


6.39(+3) 


6.42(+3) 


6.34(+3) 


6.34(+3) 


6.34(+3) 


6.37(+3) 


C0 2 


90 


8.95(+l) 


2.29(-7) 


1.24(+2) 


1.24(+2) 


1.24(+2) 


8.97(+l) 



5. Aside from O2, CO, H2CO and CO2, all species show an ex- 
cellent match to the Monte Carlo and/or master equation results. 
The worst, CO2, is only a factor of ~ 1.4 greater than the master- 
equation value; a great improvement over the rate-equations. 

Table 5 shows results for the intermediate density case. Here, 
the simple competition scheme produces a good match with all 
species except CO2, which is inaccurate by a factor similar to 
that of the low-density case. The quality of the match for all 
other species is not in general as good as for the low-density 
case, but is still quite acceptable. 

Table 6 shows results for the high density case. The results 
of the simple competition scheme are generally within an order 
of magnitude of the Monte Carlo and/or master-equation results, 
and some are rather closer; O2 is very accurate. 

In all regimes, method A provides at least a reasonable 
match to the exact results, and great improvements over the rate- 
equation method. However, there are still discrepancies that can- 
not be explained merely as computational inaccuracies between 
the different methods. Further refinements to the treatment of 
surface processes are required to match the exact results. 

6. 1 . Full, two-particle competition 

Section 5 outlined a simple competition scheme that introduced 
an efficiency to the surface production rates. The scheme as- 
sumes that if one or other particle evaporates before reaction can 
occur, then the production process is over. However, this is an 
over-simplification; the final rate must take into account all op- 
portunities for an individual, waiting particle to react before it is 
otherwise removed. Consider particle A waiting on a grain for 
its reaction partner, B. Particle B accretes onto the grain, and if 
A and B meet before evaporation of either species occurs, then 
they react. If particle A itself evaporates before reaction can oc- 
cur, then the process is over, as the part of the production rate 
associated with this scenario, Rf orm (B) ■ (N(A)), presumes the 
presence of particle A on the grain. But if particle B evaporates, 
particle A is still present and waiting for another particle B to ac- 



crete, of which there is a steady stream. If particle A can remain 
on the grain for another accretion timescale of B, without evap- 
orating, then the reaction competition process may begin again. 
Two competition processes must therefore be considered: com- 
petition between reaction and evaporation of A or B; and compe- 
tition between accretion of B and evaporation of A. Ignoring in- 
terference from other accreting particles, re-accretion may occur 
an arbitrary number of times. Thus, method A underestimates 
the reaction efficiency, particularly if reaction is slower than one 
of the evaporation rates. 

Full, two-particle competition requires two efficiency fac- 
tors, t] AB (A) and tj ab {B), depending on whether A or B is waiting 
on the grain. This gives a modified production rate: 

R m „d{AB) = Rf orm (B) ■ (N(A)) ■ t]ab(A) 

+ R form (A) ■ (N(B)) ■ r, AB (B). (19) 

To construct r] AB (A), or equally, r] AB (B), each individual compe- 
tition process is considered. Firstly, the probability of reaction 
between A and B is defined in equation (17) as C AB . 

The probability that evaporation of B takes place before re- 
action can occur, with A remaining on the grain, is: 

D AB (A) = k evap (B)/[k AB + k evap (A) + k evap (B)]. (20) 

In this case, A has a chance to react with another arriving par- 
ticle B. The probability that A remains on the grain, avoiding 
evaporation, long enough for another particle B to arrive is: 

E AB (A) = R form (B)/[R form (B) + k evap (A)]. (21) 

Hence, the probability that reaction does not occur at one oppor- 
tunity, but that particle A remains on the grain and has a chance 
to react with the next incoming particle B, is: 

F AB (A) = D AB (A) ■ E AB (A) (22) 

where > F AB (A) > 1. For a waiting particle A, the rate asso- 
ciated with the first particle B to arrive is, naturally, Rj orm {B). 
However, if reaction only successfully occurs with the second 
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Table 5. Modified rate results for intermediate density, using input values from Stantcheva et al. (2002). Values in boldface show 
agreement within 10% of the Monte Carlo or master-equation values. 



Species Stantcheva et al. Rate Eq. Modified Rates 





Monte Carlo 


Master Eq. 22211 




Method A 


Method B 


Method C 


Method D 


"Simple" 
Competition 


"Full" 
Competition 


Method B + 
Poisson Prob. 


Method C + 
H-accretion 


H 


1 


2.88(-3) 


5. 41 (-6) 


3.10(-3) 


3.10(-3) 


3.10(-3) 


2.82(-3) 


O 





1.36(-1) 


9.78(-6) 


1.10(-1) 


1.10(-1) 


1.16(-1) 


1.25(-1) 


OH 





1.35(-1) 


9.78(-6) 


1.10(-1) 


l.lO(-l) 


1.17(-1) 


1.25(-1) 


H 2 


0, 1 


7.01(-1) 


2.20(+l) 


7.55(-l) 


7.55(-l) 


7.54(-l) 


6.85(-l) 


o 2 


9.40(+3) 


9.03(+3) 


1.26(-4) 


9.36(+3) 


9.36(+3) 


9.36(+3) 


8.44(+3) 


H 2 


6.02(+4) 


5.93(+4) 


8.55(+4) 


5.77(+4) 


5.77(+4) 


5.77(+4) 


6.18(+4) 


CO 


1 


7.76(-l) 


4.15(+2) 


7.24(-l) 


7.24(-l) 


7.24(-l) 


7.97(-l) 


HCO 





1.14(-1) 


7.39(-6) 


1.06(-1) 


1.06(-1) 


1.12(-1) 


1.17(-1) 


H 2 CO 


1 


7.11(-1) 


4.24(+2) 


6.35(-l) 


6.35(-l) 


6.35(-l) 


7.27(-l) 


H 3 CO 





1.22(-1) 


7.39(-6) 


1.06(-1) 


1.06(-1) 


1.12(-1) 


1.17(-1) 


CH3OH 


5.79(+4) 


5.81(+4) 


6.38(+4) 


5.56(+4) 


5.56(+4) 


5.56(+4) 


5.78(+4) 


co 2 


6.60(+3) 


6.64(+3) 


9.50(-5) 


9.06(+3) 


9.06(+3) 


9.06(+3) 


6.84(+3) 



particle B to arrive, then particle A should have waited a fur- 
ther accretion/formation period, tf orm (B). Hence, the arrival rate 
associated with reaction of A with the second particle B to ar- 
rive is Rf orm (B)/2. Similarly, for the third B particle, the rate is 
Rf orm (B)/3, and so on. 



Using the arrival rates and efficiencies defined above, the 
overall production rate in the case of particle A waiting on the 
grain is constructed thus: 



Rmod, A {AB) = R form (B) ■ (N(A)) ■ Cab 

Rform(B) 



2 

Rform(B) 



(N(A)) ■ C AB ■ F AB (A) 
(N(A)) ■ Cab ■ F 2 AB (A) 



= R form (B) ■ (N(A)) ■ Cab Y 1 -^-t- 

n=0 

= R form (B) ■ (N(A)) ■ Cab f -^77- In [1 - F AB (A)] 
{ F AB (A) 

The efficiency is therefore defined as: 

^b(A) = - In [1 - F AB {A)] (23) 

Fab(A) 

where > 7] AB (A) > 1. As F ab (A) 0, i] AB (A) C AB - 
Calculation of t] A b(A) and t]ab(B) is trivial, and has a minimal 
effect on the run-time of the program. 

Tables 4-6 show, under the heading "Method B", the re- 
sults of the application of these efficiencies to all reactions in 
the methanol system. In fact, the modification has a discernible 
effect only on abundances of CO and H2CO, and then only in 
the low-density regime. Reactions involving these two species 
are hindered by activation energy barriers, making them uncom- 
petitive compared to hydrogen evaporation; hence, their efficien- 
cies are especially strongly affected by the change. Indeed, the 
abundances of CO and H2CO for low density, using the modi- 
fied rates, are now an exact match to the master-equation results, 
within computational accuracy. In the case of the intermediate- 
density regime, the abundances of CO and H2CO are great 
enough that the modified rates exceed the standard determin- 
istic rates, so equation (13) comes into effect. Hence, the CO 
and H2CO rates are unaffected, whichever competition scheme is 
employed. In the high-density regime, (N(CO)), (Af(H 2 CO)) » 
1, putting them in the deterministic regime. 
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Table 6. Modified rate results for high density, using input values from Stantcheva et al. (2002). Values in boldface show agreement 
within 10% of the Monte Carlo or master-equation values. 



Species Stantcheva et al. Rate Eq. Modified Rates 





Monte Carlo 


Master Eq. 23311 




Method A 


Method B 


Method C 


Method D 


"Simple" 
Competition 


"Full" 
Competition 


Method B + 
Poisson Prob. 


Method C + 
H-accretion 


H 





8.29(-9) 


3.39(-10) 


3.26(-9) 


3.10(-9) 


8.98(-9) 


7.37(-9) 


O 


1 


5.76(-l) 


4.52(-l) 


4.51(-1) 


4.5K-1) 


4.50(-l) 


4.97(-l) 


OH 


1 


5.97(-l) 


4.52(-l) 


4.51(-1) 


4.5K-1) 


4.50(-l) 


4.97(-l) 


H 2 





1.89(-6) 


8.63(-8) 


7.52(-7) 


7.15(-7) 


2.07(-6) 


1.70(-6) 


o 2 


2.81(+5) 


2.68(+5) 


2.73(+5) 


2.72(+5) 


2.72(+5) 


2.71(+5) 


2.68(+5) 


H 2 


1.79(+5) 


1.71(+5) 


2.49(+5) 


2.25(+5) 


2.26(+5) 


1.81(+5) 


1.96(+5) 


CO 


5.28(+5) 


5.23(+5) 


6.43(+5) 


5.93(+5) 


5.94(+5) 


5.14(+5) 


5.34(+5) 


HCO 





1.53(-1) 


5.44(-3) 


5.07(-2) 


4.82(-2) 


1.24(-1) 


9.93(-2) 


H 2 CO 


5. 01 (+4) 


5.12(+4) 


1.50(+3) 


2.21(+4) 


2.16(+4) 


4.63(+4) 


3.88(+4) 


H 3 CO 





3.62(-2) 


2.61(-5) 


4.10(-3) 


3.79(-3) 


2.37(-2) 


1.63(-2) 


CH3OH 


1.10(+4) 


1.17(+4) 


4.81(+0) 


1.89(+3) 


1.82(+3) 


1.16(+4) 


7.93(+3) 


co 2 


5.82(+4) 


6.01(+4) 


1.65(+3) 


2.91(+4) 


2.83(+4) 


7.40(+4) 


6.53(+4) 



6.2. Estimation of probabilities 

The modification method developed so far shows excellent 
agreement with the exact solutions, in the low-density regime; 
only O2 and CO2 show deviations greater than may be explained 
purely by computational rounding errors or minor differences in 
input values. Results for the intermediate-density regime are also 
very close, aside from O2 and CO2. Some populations that are 
less than unity are marginally different from the master-equation 
results; however, those cannot be regarded as wholly reliable, 
due to the use of cut-offs in the master-equation scheme. Other 
species with significant abundances are acceptably close to ei- 
ther the master-equation or Monte Carlo results. 

In the high-density regime, the match is less acceptable; H2, 
H2O, H2CO and CH3OH all show significant deviations from 
the exact results. What causes these discrepancies? 

In the intermediate- and high-density regimes, the average 
populations of some stochastic species get close to unity, even 
with the exact models. In the high-density case, O and OH aver- 
age populations are around 0.5. This means that the probability 
approximation of equation (10) may no longer be accurate. 

The probability that 1 or more atoms/molecules of species i 
should be present on the grains at any arbitrary moment is: 

CO 

P(i) = Y J P N (i) = l-P (i) (24) 

N=l 

where Pn(i) is the probability of population state N(i). As (N(i)) 
approaches unity, the validity of equation (10) becomes ques- 
tionable, as the probabilities of populations greater than 1 may 
become non-negligible. These probabilities are dependent on the 
degree of coupling between different population states. 

To gauge the importance of the approximation, two ex- 
treme cases may be identified in the simple hydrogen-producing 
system in which the average abundance of atomic hydrogen 
would be less than 1. Firstly, consider a case in which k HyH » 
^ace(H) >> k evap (H), i.e. H-accretion is much faster than H- 
evaporation, and the rate of the reaction H + H— > H2 is much 



faster than accretion. Here, evaporation is unimportant, as an- 
other H would accrete - and reaction occur - before evapora- 
tion should take place. The population state N(H) = 2 is short- 
lived, and so probabilites fV>i(H) are negligible. In this case, 
the approximation P(i) = (N(i)) is valid in the entire range 
> (N(H)) > 1; although, in fact, the system reaches a steady- 
state value of (iV(H)) = 0.5. The validity of equation (10) is the 
result of the strong coupling between the N = 2 and N = states, 
engendered by the fast reaction between pairs of H-atoms. 

Consider a second case, in which k evap (H) » k H ^ H » 
Racc(H). Here, evaporation dominates the destruction paths in all 
population states; hence, only population states that are adjacent 
are (strongly) coupled. Such a system may be understood as a 
queuing process of type M/M/oo (in Kendall notation); the pop- 
ulation probabilities are therefore well described by a Poisson 
distribution. In this case, the probability of finding more than 
one particle of species i at any instant is: 

P(0 = l-exp[-<AT(0>]. (25) 

Thus, there are situations in this simple system where equation 
(10) is still valid, even when (N(i)) approaches unity; but there 
are others in which it is not - equations (10) and (25) represent 
the extreme cases. It may be simplistically argued that, even for 
a large reaction network, P(i) for any species with (N(i)) < 1 
should lie somewhere between these two extremes. The precise 
value of P(i) would be dependent on the coupling between pop- 
ulation states of the entire reaction network. 

Since the explicit evaluation of discreet population-state 
probabilities requires a master-equation approach, making an ap- 
proximation to P(i) is unavoidable. The simplest solution is to 
adopt one of equations (10) and (25). For a value (N(i)) = 1, 
these extreme cases differ in P(i) by a factor of 1.58. In com- 
parison to the errors inherent in chemical models, this difference 
is probably unimportant. However, one may argue that equation 
(25) is the more generally valid: In a complex network, reactions 
between many different species may occur, so it is generally less 
likely that destruction processes be dominated by reactions be- 
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tween like species, e.g. H + H — > H2. Reaction between het- 
erogeneous particles would be more likely, so the potential for 
strong coupling between non-adjacent population states (in par- 
ticular, N(i) = 0, 2) should be mitigated. 

The incorporation of equation (25) into the model is la- 
belled "Method C". Tables 4 and 5 show that in the low- 
and intermediate-density cases, the change has only marginal, 
though beneficial, effects. However, the high-density results of 
Table 6 are much improved. The abundances of many important 
species, such as CO, CO2 and H2CO, show a very good level 
of agreement, whilst H2O and CH3OH are a perfect match to the 
exact results, within computational errors. The match for species 
with abundances less than 1 is also much improved, although not 
perfect. For high density, the results for all species are within 
35% of the exact results, and the majority are within 10%. 

6.3. Hydrogen Accretion Competition 

Although method C produces an acceptable match to the exact 
methods in each density regime, certain species, most notably 
CO2 and O2, are not so well reproduced, particularly at low and 
intermediate densities. In fact, one further competition process 
must be considered to improve these results. 

The formation of CO2 and O2 depend on the mobility of the 
oxygen atom on the grain surface. In comparison to atomic hy- 
drogen, oxygen is very slow to diffuse between binding sites, 
due to its greater diffusion barrier, and its being much more 
massive than H, making tunnelling ineffective (although efficient 
tunnelling of hydrogen is also questionable, see Section 6). The 
consideration of competition up to now has concentrated solely 
on the case of two reactive species on a grain at any one time. 
However, the reaction rate for O-dependent reactions is so slow, 
at ~ 4 x 1CT 5 s _1 (using values from Table 3), that a hydrogen 
atom may accrete before reaction occurs; it may then react with 
one or other of the reactants considered in the oxygen reaction. 

To test the influence of this competition process, terms equal 
to the accretion rate of hydrogen are inserted into equations 
(17), (20) and (21), for all reactions involving atomic oxygen. 
Reaction of atomic hydrogen with oxygen (or the other reactant, 
where applicable) is assumed to be instantaneous, if H-accretion 
takes place before the oxygen-dependent reaction can occur. 

Tables 4-6 show the results of this approach, labelled 
"Method D". In the low-density regime, the reproduction of the 
exact results is now perfect, within computational accuracy. In 
the intermediate density regime, the abundances of a number of 
species are now an exact match, or very close, e.g. H, H2, H2O, 
CO, HCO, H 2 CO, CH3O and CH3OH, whilst 2 is slightly less 
accurate than with method C. Even the worst match, CO2, is 
much improved, falling within 10% of the exact result. 

But in the high-density case, whilst some species such as CO 
show improved accuracy, many others diverge from the exact re- 
sults, including formaldehyde (H2CO) and methanol (CH3OH). 
The level of agreement is generally worse than that achieved 
with method C. This may be explained by the fact that in this 
density regime, the accretion rate of oxygen itself is very fast. 
The accretion of hydrogen atoms interferes with certain reac- 
tions, but the further accretion of oxygen acts to mitigate this 
effect, as it may itself react with the newly accreted hydrogen. 

It becomes clear that in order to devise a generalised system 
to deal with accretion-competition effects, it would be necessary 
to consider not only the accretion of every reactive species, but 
also the probability that each might react with any other accret- 
ing species, in the context of a reaction between two entirely dif- 
ferent reactants. Such a scheme could certainly be devised, but 



Table 7. Equations employed for each modification method 



Modification scheme Required Equations Notes 



Basic (H 2 ) 


10- 


13 






See Section 3.1 


Basic (continuous) 


10- 


15 








Method A 


10- 


18 








Method B 


10, 13-15, 


17- 


23 




Method C 


13- 


15, 17 


-23, 


25 




Method D 


13- 


15, 17 


-23, 


25 


See Section 6.3 



it would be extremely complex, both to implement and to fully 
understand. It would also, in all likelihood, be far more compu- 
tationally expensive than, for example, Method C. Under those 
circumstances, the value of choosing such an approach over ex- 
act methods like the master equation would be questionable. 

6.4. Slow hydrogen regime 

In addition to the parameter set shown in Table 3, 
|Stantcheva et al. (20 02) implemented the "slow" (M2) surface 
rates of Ruffle & Herbst (2000). These are defined by diffusion 
barriers in line with the values suggested by |Katz et al. (1999)| 
for the H + H — > H2 reaction. Implicit in these rates is the as- 
sumption that hydrogen-atom tunnelling through the diffusion 
barrier is inefficient; they represent a purely thermal hopping 
mechanism. Using these values, rate equations produce a per- 
fectly acceptable match to the results of the Monte Carlo tech- 
nique, in all regimes. Those same results are reproduced by even 
the simplest of the modified-rate methods presented here. The 
slower rates ensure either that the rate limit of equation (13) is 
reached, or that populations are high enough to place each reac- 
tion in the deterministic limit. This would not necessarily be the 
case in regimes with higher temperatures or smaller grains. 

7. Discussion 

It is clear that even a basic rate-modification scheme, such as ap- 
plied here to the hydrogen system of Barzel & Biha m (2007b)| 
is capable of achieving accurate results. However, application 
to the more complex water- and methanol-producing systems of 
Barzel & Biham demonstrates the need to consider competition 
between reaction and evaporation processes. In fact, the influ- 
ence of such competition is strongly dependent on the choice of 
surface parameters. The barriers to diffusion chosen by Barzel 
& Biham are unusually high, around 85 - 90 % of the barrier 
against evaporation (see Table 2); this makes evaporation espe- 
cially competitive with reaction. The lower barriers employed 
for the hydrogen-only system (Eb'.Eo — 0.69, see Table 2) allow 
little scope for evaporation to compete; see Figure 1. It is for this 
reason that the basic modification scheme is so accurate for the 
H2 system. Systems with yet lower E/, : Ed ratios would also be 
very well reproduced with the basic method. 

Indeed, the explicit treatment of evaporation-reaction com- 
petition in the methanol system of Stantcheva et al. (2002) makes 
only a small difference to the results, as evaporation is rela- 
tively slow. The greatest challenge in reproducing the low- and 
intermediate-density results presented by those authors is the 
treatment of reactions with activation-energy barriers. For such 
reactions, the competition comes not from evaporation but from 
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the accretion of other species with which the reactants may pref- 
erentially react. This situation may be effectively treated by in- 
cluding the relevant accretion terms in the competition efficien- 
cies; but, the approach is not universally applicable, as the high- 
density results demonstrate. However, method D is less neces- 
sary in that case, as method C already produces an acceptable 
match. In all regimes, method C still reproduces most abun- 
dances within 10% of the exact results, with the remainder of 
species within a few tens of percent of their true values. 

Whilst the neglect of accretion competition appears to make 
only a small difference to the results of the particular systems 
modelled here, would a change in parameters (e.g. accretion 
fluxes, diffusion barriers) lend a greater influence to such effects? 
Consider the formation of CO2 and O2, the only two species not 
accurately reproduced at low density by methods B and C. In 
that regime, formation of CO2 occurs mainly by the reaction O + 
HCO -» CO2 + H, due to the high activation energy required for 
the alternative route. Also, the CO population is small, meaning 
that HCO is formed stochastically (see Section 5.1). Hence, the 
production rate of CO2 should be R acc (0) ■ (A^(HCO)); however, 
according to equation (13), the production rate may not exceed 
the deterministic maximum of k ,Hco • {N(HCO)) • (N(O)). The 
destruction of atomic oxygen is dominated by (stochastic) reac- 
tion with atomic hydrogen, so its population is given approxi- 
mately by (N(0)} = R acc (0)/R acc (U). As a result, the forma- 
tion of CO2 falls into the deterministic regime when R acc (H) > 
ko,Hco- The implication for accretion competition is clear: hy- 
drogen accretion cannot become more competitive than reaction, 
before the deterministic limit is reached and rate equations take 
over. The competition efficiency factor will thus never fall be- 
low Co,hco = 1/2. A similar analysis yields Co,o = 2/3 for O2 
production. A parameter search that varies fluxes and diffusion 
barriers around their low-density values bears out these analy- 
ses: CO2 and O2 populations using method C vary from those 
of method D by no more than the stated factors, whilst popula- 
tions of other species are an exact match. So far as the systems 
considered here are representative, method D may be considered 
unnecessary for the attainment of reasonably accurate results. 

Some small inaccuracies are evident using method C in the 
high-density regime. A possible cause is a break-down in the 
two-particle assumption. For example, the average populations 
of O and OH in the high-density regime are both close to 0.5. 
Whilst equation (25) takes crude account of the probability of 
more than 1 of the same species being present on a grain at the 
same time, no account is taken of O and OH both being present. 
Such a possibility should lead to competition to react with, say, 
an accreting hydrogen atom. 

The Poisson probability of equation (25) is a simple ap- 
proximation. It may be that the probability distibutions of, say, 
O and OH are skewed away from the N - 1 state if an ex- 
act method is used. However, the full results of the hybrid 
master-equation technique employed by Stantcheva et al. (not 
shown here), indicate that the populations of stochastic species 
in the intermediate- and high-density regimes are very depen- 
dent on the cut-offs employed. The use of much larger cut-offs 
in their models could give populations in closer agreement with 
modified-rate method C. Method C does in fact allow for in- 
stantaneous populations, N(i), of more than 1, through its as- 
sumption of a Poisson probability distribution - equivalent to an 
infinitely large cut-off. Its main failure is in treating the more 
complex interactions between these heterogeneous populations. 
These failures are limited by the switch-over to deterministic 
rates when population states greater than 1 become probable. 



Hence, in the intermediate- and high-density regimes, the 
only species that may be reliably compared between methods are 
those with large populations, which are treated accurately by the 
Monte Carlo method. With the exception of CO2 (whose popula- 
tion discrepancies are explained above), all abundant species in 
the high-density regime are reproduced by method C to within 
8% of the Monte Carlo model values. 

8. Conclusions 

The production-rate modification schemes presented here sub- 
stitute the deterministic production rates with so-called "small- 
grain" rates, applicable when populations are very small and 
reaction rates are fast. These basic rates are easily formu- 
lated, and have been more rigorously derived elsewhere using 
a master-equation approach (e.g. Lipshtat et al. 2004). A simple 
switching/rate-limiting system is employed for each reaction, al- 
lowing the deterministic rates to be used when appropriate. Even 
this simple scheme provides a highly accurate reproduction of 
the hydrogen system examined by Barzel & Biham (2007a,b). 

The consideration of competition between surface processes 
allows the water- and methanol-producing systems of Barzel & 
Biham also to be very accurately reproduced. A further refine- 
ment is made to better estimate the probability of the presence of 
a reactant on the grains, improving accuracy in regimes in which 
average populations approach a value of 1 . Excellent agreement 
with the results of the Monte Carlo models of Stantcheva et al. 
(2002) is achieved with this scheme, labelled method C. 

The modification method is useful in that it trades off a small 
reduction in accuracy for a scheme that does not require the ex- 
plicit treatment of population states. This allows the methods to 
be easily implemented in gas-grain chemical models. The em- 
phasis on competition processes also renders the systems under- 
standable in terms of physical processes, rather than probability 
distributions or other more abstract quantities. 

The greatest practical obstacle to achieving a perfect re- 
production of the results of exact methods is the existence of 
activation-energy barriers for certain key reactions. However, the 
method C results never fall more than a few tens of percent away 
from the true populations of species produced by these routes. 
Critically, the inaccuracies necessarily fall in a parameter space 
close to the stochastic-deterministic threshold, and are therefore 
limited by the switch-over to the standard rate-equations. Even 
these inaccuracies may be overcome by judicious consideration 
of competition processes; however, a comprehensive approach 
of this sort would be difficult, and computationally expensive. 

Perfect accuracy over all parameter space can of course only 
be achieved with exact techniques such as the Monte Carlo or 
master-equation methods. The fundamental cause of inaccura- 
cies in the modification schemes is the break-down of the two- 
particle approximation that is applied, explicitly or implicitly, 
throughout this paper. The modification schemes, culminating in 
method C, act to diminish the effects of this break-down, thereby 
bridging more closely the stochastic and deterministic regimes. 

It should be stated again that the purpose of the methods de- 
veloped here is not to achieve perfect accuracy, but rather to 
achieve acceptable accuracy in the simulation of systems that 
cannot (currently) be treated using exact methods. In this aim, 
the modification schemes presented here appear promising. In 
many regimes, some of the modifications may not even be nec- 
cessary; competition between evaporation and reaction should 
only be important if large Eb'-E D ratios are assumed, or if ex- 
tremely small grains or extremely high temperatures are consid- 
ered. The most important modification to the basic (continuous) 
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method is that applied in method C. Even the use of just the basic 
(continuous) method would be an improvement over rate equa- 
tions, in any regime. Table 7 details the equations necessary to 
implement each method. 

Certain elements of the modification method are not rigor- 
ously tested by the systems used here, and should be further 
tested in future. This includes the switch-over to deterministic 
rates when populations become large; deterministic rates typi- 
cally become valid before populations of 1 or more are achieved. 
Also, the "full competition" scheme first implemented in method 
B shows only a limited effect here. The accuracy of most simple 
reaction networks should not suffer greatly by its omission. 

The methods of production-rate modification presented are 
simple enough that they may be implemented in significantly 
larger systems than those explored here, sufficient for a full treat- 
ment of gas-phase and grain-surface chemistry. In such networks 
it may also be necessary to treat photodissociation of surface 
species, being both a source of "formation" in equation (18), 
and a source of "destruction" for the evaluation of competition 
efficiencies in equations (17), (20) and (21). Photodissociation 
is, like evaporation, a deterministic process, rendering its imple- 
mentation in the relevant equations trivial. 
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